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The long range attractive force between two hydrophobic surfaces immersed in water is observed 
to decrease exponentially with their separation - this distance-dependence of effective force is known 
as the hydrophobic force law (HFL). We explore the microscopic origin of HFL by studying distance- 
dependent attraction between two parallel rods immersed in 2D Mercedes Benz model of water. 
This model is found to exhibit a well-defined HFL. Although the phenomenon is conventionally 
explained by density-dependent theories, we identify orientation, rather than density, as the relevant 
order parameter. The range of density variation is noticeably shorter than that of orientational 
heterogeneity. The latter is comparable to the observed distances of hydrophobic force. At large 
separation, attraction between the rods arises primarily from a destructive interference among the 
inwardly propagating oppositely oriented heterogeneity generated in water by the two rods. As 
the rods are brought closer, the interference increases leading to a decrease in heterogeneity and 
concomitant decrease in free energy of the system, giving rise to the effective attraction. We notice 
formation of hexagonal ice-like structures at the onset of attractive region which suggests that 
metastable free energy minimum may play a role in the origin of HFL. 


I. INTRODUCTION 

When water is confined between two large hydrophobic 
surfaces, one finds the emergence of a surprisingly long 
ranged effective attraction between the two hydrophobes. 
The attractive force grows exponentially as the distance 
between the surfaces is lowered. This exponential sepa¬ 
ration dependence of hydrophobic attraction has become 
widely known as the hydrophobic force law. 

The first direct measurement of the interaction between 
two cylindrically curved hydrophobic surfaces in aqueous 
electrolyte solutions was performed by Israelachvili in a 
series of landmark experiments using specially designed 
surface force measurement apparatus. [1-5] These exper¬ 
iments observed that the attractive force between the 
surfaces is “long-ranged” - being detectable even as far 
as 80 A to 100 A separation. It is neither sensitive to the 
type and concentration of the electrolyte, nor to the pH 
of the solution. [5] The hydrophobic force F measured 
in experiments can be empirically fitted to the following 
form 

s = Cexp (—I) (1) 

where R is the radius of the cylindrical curve, C is an ar¬ 
bitrary constant and f is the correlation length (or, decay 
length). Accurate estimate of the range of this hydropho¬ 
bic force went through a period of uncertainty. Differ¬ 
ent experimental force-measuring techniques and differ¬ 
ent methods of hydrophobization apparently recognized 
substantial diversity in length scales. There are reports 
of a measurable attractive force at separation as large as 
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3000 A. [6, 7] However, recent experiments have led to be¬ 
lieve that the attraction at separations greater than 200 A 
is due to a variety of system-dependent, nonhydrophobic 
(or only indirectly hydrophobic-dependent) effects. [8, 9] 
The attractive forces at separations less than 100 A repre¬ 
sents the truly hydrophobic interaction. The exponential 
decay of the attractive force between hydrophobic sur¬ 
faces (or, the hydrophobic force law), however, remains 
an integral feature of all experimental results. 

Despite a great deal of research over the last three 
decades, understanding of certain aspects of the origin 
and nature of interaction between hydrophobic surfaces 
have remained unsatisfactory. Theoretical studies have 
shown that size of the hydrophobic surface, degree of 
lrydrophobicity and temperature together determine the 
critical separation [10, 11] as well as the rate of evapora¬ 
tion of confined water. [12] The phenomenon is mainly 
attributed to the formation of a liquid-vapor interface in¬ 
duced by the surface, leading to cavitation at some critical 
intersurface separation. [13] Such cavitation is also ob¬ 
served near vapor-liquid coexistence for a Lennard-Jones 
fluid confined between hard walls. [14] This de-wetting 
induced collapse is generally described by density fields. 
In this context, Lum-Chandler-Weeks (LCW) theory [15] 
has been one of the most successful approaches. It in¬ 
volves the drying interface concept, and includes effects 
of orientation of the water molecules only implicitly in the 
Hamiltonian using prior knowledge of the radial distribu¬ 
tion function. LCW theory predicts a critical intersurface 
separation of 50 A for two parallel hard plates. 

The above picture basically describes water in terms 
of inward propagating density inhomogeneity originating 
from the two surfaces. However, orientational ordering 
induced by the surfaces is ignored. Indeed, as early as 
in 1959, Kauzmann introduced the concept of “clathrate- 
structure” of water surrounding the hydrophobic solutes 
to understand the hydrophobic effect in terms of entropic 
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stability. [16] Formation of such structures emphasizes the 
role of orientational ordering in water. Subsequent com¬ 
puter simulations have also found evidence of significant 
orientational ordering among water molecules around hy¬ 
drophobic residues of protein molecules. [17] In a recent 
study, Berry and coworkers have discussed how oppositely 
directed polarization from water dipoles can give rise to 
a long range attraction. [18] 

In view of the above, the one-order-parameter (den¬ 
sity) description of Lum et al. may not be adequate to 
describe ordering of water around a hydrophobic surface. 
Hydrophobic surfaces can induce orientational inhomo¬ 
geneity by disrupting the natural orientational correla¬ 
tion among water molecules. This orientational ordering 
in turn may alter many fundamental properties of wa¬ 
ter, including density. It is thus expected to play a key 
role in describing the metastable or unstable state, if any, 
between the two hydrophobic surfaces. Comprehensive 
understanding of the long range hydrophobic attraction 
is a challenging task, particularly because of subtle nature 
of the hydrophobic effect itself and its interplay with an 
array of other forces. 

Computational approaches have proved to be very use¬ 
ful in this regard as one can explore the microscopic de¬ 
tails and try to understand their manifestation in the 
macroscopic behavior. Accurate and detailed models evi¬ 
dently gives improved accuracy in computer simulations. 
However, simplified models often give insights that are not 
obtainable at relevant time scales from complex models. 
They usually capture the essence of underlying physics 
while being more comprehensible - providing insights and 
illuminating concepts, and they do not require big com¬ 
puter resources. 

In the present work, we employ the Mercedes Benz 
(MB) model of water to study hydrophobic force law. This 
model was introduced by Ben-Naim [19, 20], and later 
parameterized by Dill and co-workers to mimic water¬ 
like properties. [21-25] The merit of this model lies in 
the simple form of the interaction potential and the re¬ 
duced dimensionality. A schematic representation of the 
model is shown in Fig. 1 and details are given in Mate¬ 
rials and Methods. Several previous studies have demon¬ 
strated that MB model generally reproduces most of the 
properties of water including the density anomaly, the 
minimum in isothermal compressibility as a function of 
temperature, and the thermodynamic properties of non¬ 
polar solvation. This model has also been used to study 
the thermodynamic properties and structural aspects of 
confined water. [26-29] To study the hydrophobic force, 
we have introduced rigid parallel hydrophobic rods made 
of non-interacting 2D Lennard-Jones (LJ) circles in MB 
water. We have performed molecular dynamics (MD) sim¬ 
ulations in NVT ensemble to investigate the structure 
and dynamics of the system. 

We find that this model system exhibits a clear attrac¬ 
tive hydrophobic force law, with a correlation length £ 
of ~ 4 ct, where a is the diameter of the MB circle. The 
force is found to decay to 90% of the contact value at 



Figure 1 . Schematic representation of the two-dimensional 
Mercedes-Benz (MB) model of water. Two MB particles are 
shown, separated by a distance r tJ . Each particle has three 
H-bonding arm vectors: ik and j t respectively (fc, l = 1,2,3). 
The interparticle axis vector is denoted by iiij and the angles 
that the closest arm of each particle make with the axis vector 
are labeled <j>i and ifij. 

around ~ 10 ct. Our results and analysis provide a molecu¬ 
lar level explanation of the origin of long-range attractive 
forces between the two hydrophobic surfaces, prior to the 
eventual cavitation induced collapse. We introduced and 
calculated appropriate orientational order parameters. We 
find that the correlation length of density heterogeneity is 
significantly shorter than that of the bond orientational 
order parameter, while the latter is comparable to the 
correlation length of the hydrophobic force. Destructive 
interference between orientational heterogeneities propa¬ 
gating inwards from the two surfaces lead to a lowering 
of imposed order, and hence lowering of free energy. As 
this destructive interference increases with decreasing sep¬ 
aration between the two rods, there appears an effective 
attraction between the two surfaces. 


II. HYDROPHOBIC FORCE LAW IN MB 
WATER 

The hydrophobic force law has been inspected in com¬ 
puter simulations studies of 3D water models, but the 
validity of the same in reduced dimension, such as in the 
MB model, has not been looked at. Recently, Dill and co¬ 
workers have studied the behavior of confined MB water in 
enclosed cavities using NVT and jiVT Monte Carlo sim¬ 
ulation. [26] The force between the walls was estimated 
from the oscillating density profiles and was found to be 
attractive. The exponential distance-dependence and the 
correlation length were, however, not examined. In this 
work we calculate the pressure P cav in the confined region 
between the two hydrophobic rods by using the virial 
expression, 
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Figure 2. Pressure P on each individual hydrophobic LJ rods 
suspended in MB water increases exponentially as the distance 
between the rods ( d ) is increased. Simulation results are shown 
for different rod lengths, given by the number of LJ particles 
(rirod) constituting the rod. The solid lines are the exponential 
fittings [Eq. 5] with global correlation length £ = 3.8 ±0.1. 


where N is the number of atoms in the system, k B is 
the Boltzmann constant, T is the temperature, D is the 
dimensionality of the system, V is the system volume, r t 
is the position vector of i th particle, and f \ is the total 
internal force (arising from the interatomic interactions) 
on i th particle. The second term in the above Eq. 2 is the 
virial. 

The pressure is found to increase exponentially with the 
distance of separation between the two rods. We obtain 
the pressure at infinite separation by fitting it to the 
equation, 

-Pcav = Poo + A exp (3) 

where d is the inter-rod distance, A and B are fitting 
parameters. The effective pressure on the rods is subse¬ 
quently obtained as 

P = P cav - P x (4) 

In Fig. 2, we show the effective pressure on the rods as a 
function of inter-rod separation distance d. The length of 
the rod is given by the number of LJ particles (n rod ) con¬ 
stituting the rod. The effective pressure shows negligible 
dependence on the rod length, and decays exponentially 
with inter-rod distance. This is consistent with the exper¬ 
iments of Israelachvili et al. [2, 5] and can be fitted to the 
equation 



Table I. Fitted values of the correlation length and prefactor 
of the hydrophobic force law (see Eq. 5) observed between LJ 
rods in a system of MB particles 


n rod 

C 


5 

-0.21 ±0.01 


10 

-0.25 ±0.01 

3.8 ±0.1 

15 

-0.26 ±0.01 


20 

-0.27 ±0.01 



which is analogous to the hydrophobic force law Eq. 1. 
The solid lines in Fig. 2 show the global fitting with con¬ 
stant £. The correlation length is found to be £ = 3.8±0.1. 
The values of the prefactor C (obtained from the fitting) 
for the different rod lengths are given in Table I. The net 
force on the rods can be obtained by multiplying the effec¬ 
tive pressure with the respective length of the rods. Thus 
the correlation length £ for the effective force would not 
depend on the rod length, though the prefactors would 
change. 

The negative pressure indicates effective attraction be¬ 
tween the rods, induced by the intervening water. As 
noted by the dashed lines in Fig. 2, the pressure in the 
confined region reaches 90 % of the bulk at d « 10. Thus, 
there exists a tangible finite force of attraction between 
the two rods even when they are separated by 10 water 
diameters. 

Since the effect of the rod length on hydrophobic force 
appears negligible, we use a particular size of the hy¬ 
drophobic rod, n rod = 15 for further analysis. 

III. MICROSCOPIC ORIGIN OF THE 
HYDROPHOBIC FORCE: ORIENTATION AND 
DENSITY MAPS 

To understand the origin of long range hydrophobic 
force, we seek to visualize the underlying control param¬ 
eters - position-dependent density and orientation - of 
MB water in the region confined by the two parallel hy¬ 
drophobic rods. 

We map the orientation of MB particles at different 
inter-rod separations in the left panel of Fig. 3 and com¬ 
pare with the corresponding density maps in the right 
panel. The angles are measured from the 2 -axis, be. the 
axis normal to the rods. The angle zero represents an arm 
pointing away from the left rod. Conversely, the angle 
180° represents an arm pointing away from the right rod. 
The MB particles (with 3 arms) near the wall prefer a par¬ 
ticular orientation, with one of the arms pointing towards 
the wall. Thus, particles near the left wall prefers 60°, 
180° and 300° angles, while particles near the right wall 
prefers 0°, 120° and 240°. It means that the particles near 
the left wall has a dangling, unsatisfied H-bond pointing 
towards the wall, and vice versa. The preference is due to 


P = Cexp 
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Figure 3. Orientation maps (left panel) and density maps 
(right panel) of MB water in the region confined between the 
two hydrophobic rods, at different inter-rod separation dis¬ 
tance (d). Here, 6 is the angle made by the arms of the MB 
particles with the 2 -axis, i.e. the axis normal to the rods. The 
left rod is at 2 = 0. The y-axis is along the length of the rods, 
with the bottom of the rods at y = 0. The orientation maps 
show the probability distribution of 9 along 2 . The density 
maps show the distribution of density in the y — z plane. 


the H-bond energy gained by the particles from forming 
“good” H-bonds with their remaining two arms. Interest¬ 
ingly, this constraint induces an orientational pattern in 
the system. Oppositely oriented patterns are created by 
the two hydrophobic walls. The inward propagation of 
the opposing orientational heterogeneity causes an inter¬ 
ference as the opposite “waves” meet near the center of 
the confined region. This results in a frustration of ori¬ 
entation near the center, and the orientational pattern 
survives for large inter-rod distances - even when the 
rods are separated by 10 water diameters. 

The density maps of MB water in the enclosed regions 
are plotted alongside in the right panel of Fig. 3 for com¬ 
parison. Uniformity along the y-axis ensure that the ori¬ 
entational patterns of the MB particles are in dynamical 
equilibrium. Visual inspection shows that oscillations in 
density at the center of the inter-rod region becomes bulk¬ 
like at d = 8.0, while the induced orientational pattern 
exists even at d = 10.0 and beyond. The emergence of 
such exotic patterns both in orientation and density in 
such simple systems is remarkable and fascinating. They 
provide a new insight to the origin of long range hydropho¬ 
bic force. 


IV. DISTANCE DEPENDENCE OF DENSITY 
AND ORIENTATIONAL INHOMOGENEITY: 

QUANTITATIVE ANALYSIS 

One primarily employs density-dependent descriptions, 
such as the LCW theory to understand hydrophobic force 
law. However, as noted in the introduction, this one-order- 
parameter description might not be enough to describe a 
complex liquid like water, especially due to the decrease in 
stability of the liquid confined within the hydrophobic sur¬ 
faces. Orientational inhomogeneity may propagate longer 
distances than density variation, as we have visualised 
in the preceding section. Here, we obtain and compare 
the range of density and orientational inhomogeneity in 
terms of correlation lengths of the corresponding order 
parameters. 


A. Length scale and amplitude of density 
heterogeneity 


Density profiles of MB particles within an enclosed 
cavity were studied by Dill et al. [26]. In the present 
work, the particles located in between the rods are in 
dynamic equilibrium with the bulk. Here, we calculate 
the distance-dependent density within the cavity using 
the same standard method, 


PcavO) 


N 

Lh 


( 6 ) 


where N is the average number of particles in an interval 
of width h at distance 2 from the left rod and L is the 
length of the rods [L = (n rod — 1 )<t lj ]. The density profile 
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Figure 5. Plot of density order parameter p C a V /pbuik as a 
function of distance d between the two rods. Here, p cav is the 
average density of MB water in the confined region, while Pbuik 
is the average bulk density of MB water. The solid squares 
are the simulation results. The solid line is the fitting to Eq. 7 
to obtain the correlation length, £ = 2.4. 


the equation 


= 1 + A exp 

Pbuik \ C / 


( 7 ) 


Figure 4. Density profile p ca v(z)/pbuik of MB water between 
the two parallel hydrophobic rods at different inter-rod sep¬ 
aration d. The oscillating density along the 2 -direction show 
that the water molecules are structured into multiple layers 
between the two rods. Analogous results were obtained by Dill 
et al. [26] in a nanotube using [iVT simulation. 


within the cavity is normalized by the density of the bulk 
(Pb„ik). In Fig. 4, we show the computed profiles of water 
density [p ca v(z)/p b uik] between the two parallel rods, at dif¬ 
ferent inter-rod separations d. The profiles are consistent 
with the earlier work of Dill et al. [26] The results show 
that the water molecules are structured into multiple lay¬ 
ers between the two rods, and hence the density oscillates 
along the z-direction. We note that similar results were 
also obtained for three dimensional water models, studied 
by Choudhury and Pettitt. [10] 

The density p cav of MB water in the region enclosed by 
the two parallel hydrophobic rods gives the scalar density 
order parameter. It is obtained by time averaging over 
the simulation trajectory, and is normalized by the bulk 
density p bulk of MB water. In Fig. 5 we show the variation 
of p cav /Pbuik as a function of the inter-rod separation d, 
while each rod consists of 15 LJ particles (n rod = 15). The 
results from simulation are shown by solid squares. The 
variation is found to be exponential, and can be fitted to 


where £ is the correlation length of the scalar density 
order parameter. The fitting is shown by the solid line 
in Fig. 5. We obtain A = —1.86 and £ = 2.4. Thus, the 
correlation length for density is less than the correlation 
length of the force between the two rods. 

A useful measure to quantify the effect of the walls on 
the water structure is to find the separation at which the 
density reaches 90 % of the bulk density. This is achieved 
at d ss 8.0 in the present case, as shown by the dashed 
lines in Fig. 5. Thus, the density of MB water in the con¬ 
fined region saturates to bulk density even though there 
persists a significant finite force of attraction between 
the rods. While the exponential decay of density may be 
correlated to the exponential decay of the force at short 
length scales, the origin of the attractive force beyond 
a certain length scale (where cavity density saturates to 
bulk density) cannot be correlated with the density order 
parameter. 


B. Length scale and amplitude of orientational 
heterogeneity 

In order to quantify the extent of orientational order 
present in the confined region, one must define an orien¬ 
tational order parameter. As we have already noted, ori¬ 
entations of the MB particles are opposite to each other 
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near the two walls. Hence they induce opposing patterns, 
which propagate towards the center. Towards this goal, we 
introduce the orientational order parameter Q = (cos 30) 
to quantify the orientational order propagation. Here 0 is 
the angle made by the arms of the MB particles, relative 
to the 2 -axis, and (...) denotes the averaging. The value 
of (cos 30) can range from — 1, when all particles are per¬ 
fectly ordered with arms pointing away from the left wall, 
to +1 when all particles are perfectly ordered with arms 
pointing away from the right wall. 

The order parameter (cos 3 6) has been calculated for 
MB water within small intervals of width h at a distance 2 
from the left rod. The distance-dependent propagation of 
Q is shown in Fig. 6 at different inter-rod separations. We 
see interference of the orientational order induced by the 
two rods. Strong oscillations are observed in the central 
region. Adjacent positive and negative values indicate 
reversal of orientations, due to interference among the 
inwardly propagating oppositely oriented heterogeneity. 
It is interesting to note that fluctuations in Q persist for 
longer distance than that in density profile. 


C. Decay of global orientational order 

A global order parameter within the cavity could be 
helpful for comparison with the scalar density order pa¬ 
rameter and formulation of a phenomenological theory. 
Here we use the six-fold bond orientational order 06, 
which is defined as 
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( 8 ) 


where N n b is the number of nearest neighbors of a tagged 
particle, Okj is the angle of the bond that the j th neigh¬ 
bor makes with the tagged particle and N is the total 
number of MB particles in the confined region. We had 
earlier shown that 06 is an effective order parameter to 
distinguish the honeycomb solid and liquid state of the 
MB liquid [30]. 

In Fig. 7, we show the distributions of 06 in the confined 
region, at different inter-rod separations and compare it 
with that of bulk MB liquid. The long tails of the distri¬ 
butions reveal the ordering in the confined region. While 
bulk-like behavior is approached with increasing inter-rod 
separation, it is interesting to note the significant ordering 
in the confined region even at d = 20.0. In Fig. 8 we show 
the variation of 06 as function of inter-rod separation d. 
The results from simulation are shown by solid squares. 
Initially, when there is drying in the cavity, there are not 
enough particles for any ordering. At d = 3.5 we find 
maximum ordering in the confined region, similar to that 
of an honeycomb solid, inside the cavity. This might be 
due to a metastable state of the MB system (akin to low 
density liquid in 3D water). The decay in orientational 
order is found to be exponential and can be fitted to the 


Figure 6. Oscillations in orientational order parameter 
(cos 38) along the 2 -axis of separation of the two parallel hy¬ 
drophobic rods at different inter-rod separations d. Results 
from simulations are the points, whereas the solid lines are 
merely visual aid. 


following form 


06 - 0oo 



(9) 


where 7 is the correlation length of the orientational order. 
The fitting is shown as a solid line in Fig. 8 . We obtain, 
A = 0.5, 7 = 3.4 and 0oo = 0.04. It is interesting to 
note that the correlation length for orientational order is 
much longer ranged than the density order (0 = 2.4), and 
comparable to that of the hydrophobic force (£ = 3.8). 


D. Spatio-temporal correlation of orientational 
relaxation 


The orientational relaxation timescale provides further 
insight to the dynamics of orientation. The orientational 
relaxation of the MB particles is charcterized by the ori- 
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Figure 7. Probability distribution of six-fold bond orienta¬ 
tional order ipa of MB liquid confined between the two hy¬ 
drophobic plates. The distribution at different inter-rod sep¬ 
aration d are color coded, as shown in the legend. The long 
tail of the distributions of confined MB liquid, as compared to 
the distribution in bulk, indicate occurrence of ordered liquid 
within the cavity. 
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Figure 8. Average six-fold bond orientational order ipe of 
MB water in the region between the two parallel hydrophobic 
rods, as a function of distance d between the two rods. The 
solid squares are the simulation results. The solid line is the 
fitting to Eq. 9, showing that the correlation length 7 = 3.4 
propagates longer distance than density order parameter, and 
is comparable to that of the hydrophobic force (see Fig. 2). 


entational correlation function C(f), 
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Figure 9. Plots of average orientational relaxation time (r) 
of MB water along the 2 -axis of separation of the two parallel 
hydrophobic rods at different inter-rod separations d. Results 
from simulations are the points, whereas the solid lines are 
merely visual aid. 


where 9(t) is the angle made by the arms of the MB 
particles, relative to the z-axis at time t, N is the number 
of particles in the region considered, and (...) denotes the 
ensemble average. We calculate C(t) for the MB particles 
within small intervals of width h at a distance z from 
the left rod. We have considered time average instead of 
the ensemble average, since they are equivalent because 
of ergodicity. The orientational relaxation is found to 
depend on z. We fit it to a biexponential form 

C(t) = 01 exp + (1 — ai) exp (^~^j (H) 

where the time constants for the decay of the orientational 
correlation are given by 77 and t 2 , with amplitudes a\ and 
(1 — ai) respectively. The average relaxation time (r) is 
subsequently obtained as (r) = a\T\ + 0272 where a 2 = 
(1 — ai). The variation of (r) with distance at diffferent 
inter-rod separations d is shown Fig. 9. The orientational 
relaxation is found to be faster near the hydrophobic 














surface, while it reduces as we approach the bulk-like 
liquid near center of the confined region. For large inter 
rod separations, the average relaxation time saturates to 
the bulk value of r bu i k = 11.3 near the center. 


V. FREE ENERGY SURFACE FOR 
CONFINEMENT INDUCED TRANSITION: 

POSSIBLE ROLE OF A METASTABLE PHASE 

As the density of the liquid in the confined region de¬ 
creases below the bulk value and ice-like orientational cor¬ 
relation increases (see Fig. 6-9), additional contribution 
to the stability of these low density molecular arrange¬ 
ments may come from the presence of nearby metastable 
phases. [31, 32] In the present problem, the low density 
metastable phase is the hexagonal 2D ice (honeycomb). 
For the 3D water, nearby metastable phases are ice lh and 
the low density liquid(LDL)-like arrangement. For MB 
particles, there is an additional metastable phase which 
is oblique phase. [30, 31, 33] However, oblique phase is 
quite far off in the parameter space and we do not see 
any evidence of oblique phase. However, we see strong 
evidence of the presence of ice-like configurations (see 
Fig. 10). As discussed elsewhere [31], metastable minima 
in free energy can lower the surface energy and thereby 
play an important but “hidden” role in the stability of the 
confined liquid, and hence provide additional attraction 
at some separation. 

As mentioned earlier, existing theories of hydropho¬ 
bic force are based primarily on number density as the 
sole order parameter that varies across the system as 
we travel from one surface to another. As a result, they 
can miss important physical constraints coming explic¬ 
itly from orientation of the water molecules, such as the 
clathrate-like structures with pronounced orientational 
ordering. [16, 17, 34] It becomes particularly important 
when the separation is somewhat large, may be of the or¬ 
der of 8-10 molecular diameters wide in the present MB 
system. At such separations, the density inhomogeneity 
become rather small but orientational heterogeneity may 
still be significant, as we indeed find in our simulations. 
This suggests that orientational heterogeneity is impor¬ 
tant and needs to be taken into account. The crucial role 
of orientation was earlier pointed out by Tanaka [35], in 
case of crystallization, quasi-crystal formation, glass tran¬ 
sition, etc. in 3D liquids. In fact, for water-like systems, 
orientation could be more important of the two as density 
is slaved to orientation imposed by H-bonding. 

Physical picture becomes more interesting as we move 
to larger d due to enhanced population of 6-membered 
rings. In Fig. 10, we show some representative snapshots 
of the simulations at different inter-rod separations. We 
find an increase in the presence of hexagonal ice-like rings 
at intermediate separation distances (d = 8—12) between 
the rods. As already discussed, this corroborates with the 
argument that the system experiences the metastable free 
energy minimum due to the ice phase. Clearly, such a min- 



(a) d = 4.0 (b) d = 6.0 


(0 d = 8.0 (d) d = 10.0 

Figure 10. Representative snapshots from the simulation, 
showing the region of hydrophobic confinement, at different 
inter-rod separations d. The whole box containing 2160 MB 
particles is not shown, instead a specific region has been fo¬ 
cused for clarity. 


imum can help lowering surface free energy if the system 
permits orientational arrangement close to ice. Coopera- 
tivity of water molecules induced by the hydrophobic sur¬ 
face has been observed in 3D water models as well. [17, 34] 

In Fig. 11 we show the calculated free energy surface 
corresponding to the melting of MB system at P = 0.19 
and T = 0.15. The minimum at high p and low ipg corre¬ 
sponds to the liquid phase, whereas the honeycomb phase 
(akin to low density ice) creates another minimum at lower 
p and higher i/jq. The increase in orientational order and 
decrease in density, as observed in our simulations and 
visualized in the snaps (Fig. 10) may correspond to this 
metastable minimum. While LCW theory considers the 
cavitation and gas-liquid transition for collapse of large 
hydrophobic surfaces, this metastable minimum should be 
considered for longer ranged attraction. 

The above conclusions are substantiated by the sharp 
rise in fluctuation of density and orientational order as 
shown in Fig. 12. We show the variance of density and 
orientation order parameters at different inter-separation 
distances in Fig. 12a. The crossover at d = 4.0 suggests in¬ 
creased fluctuation and a corresponding free energy mini¬ 
mum, which might be due to metastable hexagonal phase. 
Another interesting way to explore the underlying free 
energy is to study the susceptibility for the different prop- 
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Figure 11. Free energy surface of the Mercedes Benz system 
at P = 0.19 and T = 0.15. Here, p is the density of the liquid 
and tpe is the six-fold bond orientational order. It shows the 
liquid phase at high p and low ip&, while the honeycomb phase 
(akin to low density ice) appears at lower p and higher ipe. 




d 


d 


Figure 12. Plots of order parameter fluctuations and suscep¬ 
tibilities. (a) Fluctuations in orientational order parameter 
(er^ 6 ) and cavity density (er Pcav ). (b) Compressibility (Xn oav ) 
of the MB liquid in the confined region, and comparison with 
susceptibilities of orientation (x-ipe) as we U as cavity density 
(Xpcav)■ The susceptibilities are defined in Eq. 12. 


erties, defined as, 


XA = 


(A 2 ) - (A) 2 

(A) 


(A) 


( 12 ) 


where A is the property under consideration. When we 
consider the number of MB particles in the cavity (n cav ), 
then Xn cav gives the compressibility of the system. As 
shown in Fig. 12b, the compressibility goes through a 
maximum at d = 4.0, suggesting a metastable minimum. 

A phenomenological description of hydrophobic force 
should therefore have at least two order parameters, 
namely density and orientational order (which could be 
bond-orientational order depending on the system). In 
this case of MB liquid, density between the two phases 
(liquid and honeycomb) differs by ~ 20 %, whereas the 


bond orientational order differs by ~ 75 %. A Ginzburg- 
Landau free energy functional can be written in terms 
of two order parameters, one conserved (number density) 
and one non-conserved (orientational density), to account 
for the non-local inhomogeneity, 

Ftot {V>6,p} = F-i/’e {^e} + F p {p} + F int {ipQ, p] (13) 

The first two terms on the right hand side correspond to 
each order parameter, and the third term considers the 
interaction between them. 

The free energy per particle of the MB system should 
be maximum near the walls due to higher ordering, and 
should gradually reduce towards the bulk. One has to 
consider that the ordering is induced on both sides of 
the rod, so that the excess total free energy of the sys¬ 
tem for infinitely separated rods is 4Fq“, where Fq“ is 
the maximum free energy obtained by allowing the order 
to propagate on any one side up to infinity. As the two 
rods are brought closer, the opposing orders interfere, cre¬ 
ating a frustration and reduction of order at the center 
(see Fig. 3). Thus, the total free energy of the system 
decreases up to the extreme limit of 2Fq“ when the two 
rods collapse. Asymptotically, the total change in free en¬ 
ergy in going from infinite separation to collapsed state is 
—2Fq“. This explains semi-quantitatively the attraction 
between the two rods at longer length scales. 


VI. CONCLUDING REMARKS 

Long range hydrophobic attraction plays a key role in 
several biological self-organizing processes, such as protein 
folding, assembly of membrane structures, ligand binding 
to hydrophobic patches on target proteins, micellar ag¬ 
gregation, lipid bilayer formation, etc. [36-39] It is also 
crucial to phenomena such as mineral flotation, wetting, 
coagulation, and many processes that involve surfactant 
aggregation, e.g. detergency. Deciphering the hydropho¬ 
bic interactions can also have profound implications for 
drug design and delivery to specific target proteins in a 
cell. [40] 

The present work demonstrates that a significant, at¬ 
tractive force can exist between two hydrophobic rods in 
a simple 2D mmolecular liquid. The system exhibits an 
exponential increase in the attractive force as the two 
rods are brought closer together. In addition, even the 
estimate of the range of the force is in rough agreement 
with experimentally observed range. 

Although the MB model does not possess the complex¬ 
ity of interactions experienced by real water molecules, 
it mimics many of the important anomalies and proper¬ 
ties of water. The simplicity of the present model and the 
ease of visualization can therefore be harnessed to gain im¬ 
portant insights needed to understand the puzzling (and 
often controversial) aspects of water. These are sometimes 
obscured by the complexity of the 3D models. 

Microscopic analysis of the structure and orientation of 
the system reveals several important issues. The length 
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scale of density fluctuations is found to be shorter in range 
than the length scale of hydrophobic force between the 
rods. In contrast, orientational heterogeneity propagates 
to longer lengths with a correlation length that compares 
well with that of the hydrophobic force law. Hydrophobic 
force law is found to be a consequence of the minimum 
frustration of the extended H-bond network the frustra¬ 
tion being imposed by the hydrophobic surface. Destruc¬ 
tive interference between the orientational heterogeneity 
induced by the two surfaces creates a metastable ice-like 
phase in the central region, thereby lowering the free en¬ 
ergy of the system. 

With the insight gained from the MB model, it would 
be interesting to see an extension of this work to 3D water 
models with specific attention to orientation of the wa¬ 
ter molecules, which has largely been overlooked in past 
investigations. A detailed theoretical analysis with the ori¬ 
entational order parameter may be rewarding, and could 
complement our present understanding. On a slightly dif¬ 
ferent note, pattern formation in nature has been a fasci¬ 
nating and demanding held of research and the interfer¬ 
ing orientational patterns appearing in this system (see 
Fig. 3) could be an excellent model to study. Role of ori¬ 
entation in the hydrophobic force law indicates that such 
interference as observed here may play an important role 
in directing protein association and self-assembly where 
attraction between extended hydrophobic surfaces is be¬ 
lieved to be involved. 


MATERIALS AND METHODS 

A. Mercedes Benz model 

In the Mercedes-Benz (MB) model, water molecules are rep¬ 
resented as two dimensional disks with three symmetrically 
arranged arms, separated by an angle of 120°, as in the Mer¬ 
cedes Benz logo. Each molecule interacts with other through 
an isotropic Lennard-Jones (LJ) and an anisotropic hydrogen- 
bond (HB) term in a completely separable form as 

U (Xi, X,') = (r«) + Hhb (Xi, X,) (14) 


where e LJ and ct lj are the well-depth and diameter for the 
isotropic LJ interaction. The anisotropic HB part of potential, 
I/hb is defined as 

U HB (X,.X,) = 

3 

£hb G(fij — ^hb) 'y ' G [ik ■ Uij — 1^ G (y jj ■ Uij + 1 'j 

k,l =1 

( 16 ) 

where G(x) = exp(—x 1 2 /2<r 2 ). This form of the anisotropic 
potential ensures that neighboring water molecules form ex¬ 
plicit H-bonds (have favorable energy) when the arm of one 
molecule aligns with the arm of another - the strongest bond 
occurs when they are perfectly collinear and are opposite in 
direction. The unit vector ik represents the fc th arm of the 
i th particle (fc = 1,2,3) and iiij is the unit vector joining 
the center of particle i to the center of particle j. We have 
used the optimal parameters of the MB model that are known 
to reproduce the properties of water, e H B = —1.0, r H B = 1.0, 
£lj = 0.1, cr L j = 0.7 and <t H b = 0.085. 


B. Simulation details 

In this study, we use MB model of water. The rigid hy¬ 
drophobic rods are made up of non-interacting 2D Lennard- 
Jones (LJ) particles. The LJ particles were separated by a 
distance of 0.7<t H b. For the interaction between the LJ parti¬ 
cles and MB particles, we considered the same parameters as 
used by Debenedetti et al. [12] in their 3D simulation, but 
in reduced units, Clj-mb = 0.006 and ct L j-mb = 1.2. We incor¬ 
porated the 2D MB force field in LAMMPS for carrying out 
Molecular Dynamics (MD) simulation. The forces and torques 
were calculated using similar techniques as used earlier [41] 
for the 3D model. 

We generate an initial configuration of 2160 MB particles in 
a honeycomb lattice. We perform equilibration at temperature, 
T = 0.2 and pressure P = 0.19 to obtain a liquid configura¬ 
tion. We insert two rods at required distance by replacing 
any MB particles, if within contact distance. The system was 
equilibrated for 10 s steps at constant temperature and volume 
(NVT ensemble), with each time step r = 0.007. The produc¬ 
tion run was carried out for 10 7 steps at constant temperature 
and volume. We used Nose-Hoover thermostat to keep the 
temperature bath at T = 0.2. 


where X, denotes vectors representing both the coordinates 
and orientation of the i th particle and rtj the distance between 
the centers of two particles. The notations are summarized in 
Fig. 1. The LJ term U L j is defined as 


Glj ( Tij ) 




(15) 
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